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Abstract. In this paper we provide an extensive classification of one and two dimensional dif- 
fusion processes which admit an exact solution to the Kolmogorov (and hence Black-Scholes) 
equation (in terms of hypergeometric functions) . By identifying the one-dimensional solvable pro- 
cesses with the class of integrable superpotentials introduced recently in supersymmetric quantum 
mechanics, we obtain new analytical solutions. In particular, by applying supersymmetric trans- 
formations on a known solvable diffusion process (such as the Natanzon process for which the 
solution is given by a hypergeometric function), we obtain a hierarchy of new solutions. These 
solutions are given by a sum of hypergeometric functions, generalizing the results obtained in the 
paper "Black-Scholes Goes Hypergeometric" 0. For two-dimensional processes, more precisely 
stochastic volatility models, the classification is achieved for a specific class called gauge-free 
models including the Heston model, the 3/2-model and the geometric Brownian model. We then 
present a new exact stochastic volatility model belonging to this class. 



1. Introduction 

For most mathematical models of asset dynamics, an exact solution for the corresponding 
Kolmogorov & Black-Scholes equation is usually not available: there are, however, a few notable 
exceptions. The known solutions for local volatility models arc the constant clastic of variance 
(CEV) including the classical log-normal Black-Scholes process. For the instantaneous short 
rate models, there arc the CIR process ^U] (Bessel process) and the Vacicek-HuU- White process 
[T!^ (Ornstein-Uhlenbeck process). For stochastic volatility models, the known exact solutions are 
the Heston model the 3/2-modcl and the geometric Brownian model ^HI- These analytical 
solutions can be used for calibrating a model quickly and efficiently or can serve as a benchmark 
for testing the implementation of more realistic models requiring intensive numerical computation 
(Monte-Carlo, partial differential equation). For example, the existence of a closed- form solution 
for the price of a vanilla option in the Heston model allows us to quickly calibrate the model to 
the implied volatilities observed on the market. The calibrated model can then be used to value 
path-dependent exotic options using, for example, a Monte-Carlo methodology. 

In this paper, we show how to obtain new analytic solutions to the Kolmogorov & Black-Scholes 
equation, which we refer to as KBS throughout the rest of the paper, for Id & 2d diffusion processes. 
In order to get to our classification, wc first present a general reduction method to simplify the 
multi-dimensional KBS equation. Rewriting the KBS equation as a heat kernel equation on a 
Riemannian manifold endowed with an Abclian connection, we show that this covariant equation 
can be simplified using both the group of diffeomorphisms (i.e. change of variables) and the group 
of Abelian gauge transformations. In particular for the models admitting a flat Abelian connection, 
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there always exists a gauge transformation that ehminatcs the Abehan connection of the diffusion 
operator. 

In the second part we apply the reduction method, previously presented, to a one-dimensional, 
time-homogeneous diffusion process. Modulo a change of variable, the metric becomes flat and the 
Abelian connection is an exact one-form for which a gauge transformation can always be applied. 
Using these two transformations, the resulting KBS equation becomes a Euclidean Schrodinger 
equation with a scalar potential. Extensive work has already been done to classify the set of scalar 
potentials which admit an exact solution. In particular, using a supersymmetric formulation of 
the Schrodinger equation which consists in doubling the BKS equation with another equation, 
wc show how to generate a hierarchy of new solvable diffusion processes starting from a solvable 
diffusion process (for e.g. a Natanzon potential [221 )■ In this context, the local volatility function 
is identified with a superpotential. Applying supersymmetric transformations on the Natanzon 
potential (which is the most general potential for which the Schrodinger equation can be reduced to 
either a hypergeometric or a confluent equation) , we obtain a new class of solvable one-dimensional 
diffusion processes which are characterized by six parameters. 

The classification of one-dimensional time-homogeneous solvable diffusion processes for which the 
solution to the KBS equation can be written as a hypergeometric function has been achieved in 
[21 El 01 IHl 115) using the well known Natanzon classification. The application of supersymmetric 
techniques to the classification of solvable potentials for the Schrodinger equation has been reviewed 
in [H] where a large number of references can be found. For the Kolmogorov & Fokker-Planck 
equations, one can consult |17l I15j . 

In the last part we pursue this classification for stochastic volatility models which admit a flat 
Abelian connection: we refer to these as gauge-free models. Surprisingly, this class includes all 
the well known exact stochastic volatility models (i.e. the Heston model, the 3/2-model and the 
geometric Brownian model). For this gauge-free models, we reduce the two-dimensional KBS 
equation to a Euclidean Schrodinger equation with a scalar potential. Then, we present a new 
exact stochastic volatility model which is a combination of the Heston and 3/2-models. 



In this section, we explain how to simplify the KBS equation. This reduction method will be used 
in the next section to classify the solvable one and two dimensional time-homogeneous processes. 
This method is already well known for a one-dimensional process and is presented in [3 1191 12(J| for 
example. However, the extension of this method to multi-dimensional diffusion processes requires 
the introduction of differential geometric objects such as a metric and an Abelian connection on a 
Riemannian manifold, as we will presently explain. 

Let us assume that our time-homogeneous multi-dimensional diffusion model depends on n random 
processes Xi which can be either traded assets or hidden Markov processes (such as a stochastic 
volatility a or an instantaneous short rate r). Let us denote x = {xi)i^i^... with initial conditions 
a = (ai)i=i,... ,n. These variables Xi satisfy the following stochastic differential equations (SDE) 



with the initial condition Xi{t = 0) = a^. The no-arbitrage condition implies that there exists an 
equivalent measure P such that the traded assets are (local) martingales under this measure. For 
P, the drifts 6; are consequently zero for the traded assets. Note that the measure P is not unique 



2. Reduction method 



(2.1) 
(2.2) 



dx' =b\x)dt + a\x)dWi 
dWidWj = pij{x)dt 
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as the market is not necessarily complete. Finally, the fair value of a (European) option, with 
payoff /(x*) at maturity T, is given by the discounted mean value of the payoff / conditional on 
the filtration Tt generated by the Brownian motions 

(2.3) C{a,t,T)^¥F[e-S"-^'''f\Tt] 

with rs the instantaneous short rate. This mean-value depends on the probability density p{xi, T\a, t) 
which satisfies the backward Kolmogorov equation {t = T — t, di = ^-) 

a). In this equation we have used the Einstein convention 
summed. 

can show that the fair value C for the option satisfies the 



with initial condition p(t = 0) = S{x — 
meaning that two repeated indices are 

Using the Feymann-Kac theorem, one 
Black-Scholes equation {di = 



(2.5) — = b'd,C + g'^d.fi - r{a)C 
with initial condition C(t = 0, a) = /(a)- 

In the following, the partial differential equations (PDE) 12.412.5(1 will be interpreted as the heat 
kernel on a general smooth n-dimensional manifold M (without a boundary) endowed with a metric 
Qij. The inverse of the metric g^^ is defined by 

(2.6) = \p,ja^a^ 
and the metric (p*-' inverse of py ) 



(2.7) <?y = 24- 

a a-' 

The differential operator 

(2.8) D = Vd, + g'^d^j 

which appears in 1(2.4(1 is a second order elliptic operator of Laplace type. We can then show that 
there is a unique connection V on £, a line bundle over M , and a unique smooth function Q on 
M such that 

(2.9) D = g'^V,Vj + Q 

(2.10) = g-ii^^+A^)g^g'^^J+AJ) + Q 

with g = det[gij]. We may express the connection and Q as a function of the drift 6j and the 
metric gij by identifying in 1(2.10(1 the terms di and dij with those in ((2.8(1 . We find 

(2.11) A' = i(6'_g-^a^.(gi/V^)) 

(2.12) Q = g'^Aaj - b,A - djAi) 
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Note that the Latin indices ij ■ ■ ■ can be lowered or raised using the metric gij or its inverse g''^ . 
For example Ai = gtjA^ = ^{gikb'' - \diln{g) - g.^pOkg^'^)- 

Using this connection, (|2.4|) can be rewritten in the covariant way, i.e. 

(2.13) '^^ ^ ^pi^' ^) 

If we take Ai ^ ,Q = then D becomes the Laplace-Beltrami operator (or Laplacian) A = 
g~^dig2g^^dj. For this configuration, (|2.13|) will be called the Laplacian heat kernel. Similarly, 
the Black-Scholes equation (|2.5|l can be rewritten 

(2.14) ^C{a,T) = {D-r)C{a,T) 

The Heat kernel equation can now be simplified by applying the actions of the following groups. 
t> The group of diffeomorphisms Diff{A4) which acts on the metric gij and the connection At by 

(2.15) r9^J = 9pkdJP{x)d,f{x) 

(2.16) /M, = ApdJP{x) , / e DzffiM) 

t> The group of gauge transformations Q which acts on the conditional probability (and the fair 
value C) by 

(2.17) p' = e'^(^'^)-x("^="'^=°V 

(2.18) C = e'^^^'^^C 

Then p' (C) satisfies the same equation as p (C) (|2.13l) only with 

(2.19) A', = A-d,x 

(2.20) Q' = Q + drX 

The constant phase e'^^^^"''^^"^ has been added so that p and p' satisfy the same boundary condition 
at T = 0. The above transformation is called a gauge transformation. If the connection A is an 
exact form (meaning that there exists a smooth function A such that Ai = 9,; A), then by applying 
a gauge transformation, we can eliminate the connection so that the heat kernel equation for p' 
(or C) has a connection equal to zero. It can be shown that for a simply-connected manifold, the 
statement "A is exact" is equivalent to J- ~ 0, where J- is the 2-form curvature given in a specific 
coordinate system by 

(2.21) T,, = d,Aj - d,A, 

In the following, we will restrict our classification to those processes for which = 0, meaning 
there exists a gauge transformation such that the transformed connection vanishes. The operator 
D reduces in this case to the symmetric operator _D = A + Q for which we can use an eigenvector 
expansion. 



Spectral decomposition 
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This spectral decomposition is valid only if the symmetric operator D is an (unbounded) self- 
adjoint operator or admits a self-adjoint extension. This will depend strongly on the boundary 
conditions. In order to show that D is self-adjoint or admits self-adjoint extensions, we can use 
the deficiency indices technique introduced by Von Neumann ^ f |23| and see |S] for a pedagogical 
introduction) . 

After proving that Z? is a sclf-adjoint operator or admits a self-adjoint extension, the conditional 
probability (or the fair value) can then be expanded over a complete basis of orthonormal eigen- 
vectors 4>n{x) 

(2.23) p{x, T, Xo) = ^ e^^"^(t)n{x)(t)n{xo) 

n 

with D(f)n{x) = En4>n{x). The discrete summation over n can also include a continuous summation 
according to a specific measure n{E) if the spectrum contains a continuous part. 

3. ID Time-homogeneous diffusion models 

In the next section, we apply the general reduction method, presented previously, to the one- 
dimensional KBS equation. Similar reduction to a Schrodinger equation with a scalar potential 
(without any references to differential geometry) can be found in [7| El [201 . We then find the 
supersymmetric partner of this Schrodinger equation and show how to generate new exact solutions 
(for Vanilla options). 

3.1. Reduction method. Consider a one-dimensional, time-homogeneous diffusion process with 
drift 2 



(3.1) df = fi{f)dt + <7if)dW 



Before explaining this technique, let us recall some definitions. 
An operator {H,'D{H)) defined on a Hilbert space Ti. is said to be densely defined if the subset T>{H) is dense in Ti., 
i.e. for any S H one can find in 'D(H) a sequence which converges in norm to </)■ 

The domain 'D(H^) of an adjoint operator of an (unbounded) operator H with dense domain T) is the space of 
functions V such that the linear form ip — > (i/i, Hij>) is continuous for the norm of Ti.. Hence using Riesz' theorem, 
there exists a unique <p' such that (■(/), = (tp' ,<j>) with (., .) the scalar product on the Hilbert space. By definition, 
we set ifti/) = ^' . 

An operator H is called symmetric if for all (f>,ip €. T>{H), we have {H4>,ip) = {■>p,Hxlj) . H is self-adjoint if 
additionally T>{H^) = ^{H). Let us assume that {H,'D(H)) is densely defined, symmetric and closed with adjoint 
(H^ ,'D{H^)). The deficiency indices are then defined by 

(2.22) n± = Ker{H'^ ^ ild) 

Theorem 2.0.1. For an operator with deficiency indices (n_,n+), there are three possibilities: 

(1) If n— = n_|_ = 0, then H is self-adjoint (necessary and sufficient condition) 

(2) //n_|_ = n_ = n > 1, then H has infinitely many self-adjoint extensions, parameterized by the unitary group U(n). 

(3) If n— ^ n_|_, then H has no self-adjoint extension. 

Note that if the deficiency indices are given by (n, n) then the spectrum is discrete. Moreover, if H admits a 
self-adjoint extension, the resulting conditional probability is not unique but depends on the boundary conditions 
which are parameterized by the unitary group U(n). 

■^The time-dependent process df = ii(f)A^ (t)dt -|- a{f)A(t)dW is equivalent to this process under the change of 
local time t' = A'^{s)ds 
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This process has been used as the basis for various mathematical models in finance. If / is a traded 
asset (i.e. a forward), the drift vanishes in the forward measure and we have a local volatility model 
where we assume that the volatility is only a function of /. The one-dimensional process is not 
necessarily driftless as the random variable is not a traded asset as it is the case for an instantaneous 
short rate model, or a model of stochastic volatility. 

In our framework, this process corresponds to a (one-dimensional) real curve endowed with the 
metric 5// = ^ttw- For the new coordinate 




the metric is flat: gss = 1- The Laplacc-Beltrami operator therefore becomes A= . 

Using the definition 12.11|1 . (|2.12|) . we find that the Abelian connection Af and the function Q are 
given by 



(3.3) Af^~ldM'^if))+ "^^^ 



(3.4) Q^-M{f)<r {f)--<T'Uf ) 



4V V. 2 v.. ; 2 a{f) 2a^f) 

In this case, by applying a gauge transformation on the conditional probability p, P ^ ^^^e^p 
with 

then the connection vanishes and P satisfies a heat kernel with a scalar potential Q{s) (in the s 
flat coordinate) 



(3.6) a,F(s,T) = (5,2 + Q(s))P(s,r) 

The solution P has been scaled by the (constant) factor in order to have the initial condition 
limr-,oP{s,T) ~ (5(s). Moreover, Q is given in the s coordinate by 

(3.7) Q^\ilni.))\s)-hiHa))\s)f ^'^'^ .^MsV^s) ,\s) 



where the prime ' indicates a derivative according to s. 

Example 3.1.1 (quadratic volatility process). Let us assume that / satisfles a driftless process (i.e. 
Ii{f) = 0). The Black-Scholes equation reduces to the heat kernel on M if Q{s) = constant (i.e. Q{s) 
is zero modulo a time-dependent gauge transformation) which is equivalent to cr(/) = a/^ -I- Pf + 'J 
(i.e. the quadratic volatility model, also called the hyperbohc model [T^). 



Example 3.1.2 (CEV process). For the CEV process df ~ f^dWt fJ,{f) — 0, the potential is 

/3(/3-2) 

4(1-^)2^2 



Q(«) ^ ifi^A for /3 ^ 1 and Q{s) = -i for (3^1. 
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If the Hamiltonian D = {—d^ — Q) is self-adjoint or admits self-adjoint extensions, the spectral 
decomposition can be used and the conditional probability can be decomposed using a complete 
basis of orthonormal eigenvectors (j)n{s). 

(3.8) P(s,r)=:^0„(s)0t(so)e-^"- 

n 

with (j)n{s) satisfying 

(3.9) D^^{s) = E^^^{s) 
Boundary conditions 

For a one-dimensional diffusion process, one don't need to use the deficiency indices technique as the 
complete classification of the boundary conditions is given by Feller's classification. More precisely, 
for a ID time-homogeneous diffusion process, the boundary falls into one of the four following 
types: regular, entrance, exit or natural. For entrance, exit or natural, no boundary conditions 
are needed. For a regular boundary, the conditional probability is not unique but depends on 
the boundary conditions. It corresponds to the case when the deficiency indices are (n,n). The 
boundary classification depends on the behavior of the following functions S{c, d) = s{f)df, 
M{c,d) = J^m{f)df, S(e) = lirric-^e !^ ^{c, f)m{f)df and iV(e) = lirric^e S{x,d)m{x)dx with 
s{f ) = e "t"^)^ and m{f ) = ,s , ,x (see Table 1 below). 



Boundary type 


S{t 




M{ 


e,d) 




N{e) 


Regular 


< 


oo 


< 


oo 


< oo 


< oo 


Exit 


< 


oo 




oo 


< oo 


< oo 


Entrance 




oo 


< 


oo 


< oo 


< CO 


Natural 


< 


oo 




oo 


= oo 


= oo 






oo 


< 


oo 


~ oo 


= oo 






oo 




oo 


= oo 


= oo 



Table 1. Feller's Classification. 



Example 3.1.3 (Vanilla option with constant interest rate). We have that the forward / satisfies a 
driftless process (i.e. ^{f) = 0). The value at t of a European option (with strike K and expiry 
date T) is then given by {t = T — t) 



(3.10) C{ft,K,T) = e"'-^ / (/ - K)p{f,T\ft)df 

Jk 

Doing an integration by parts, or equivalently, applying the Tanaka-Meyer formula on the payoff 
{St — K)~^ ^ we can show that the vanilla option C can be rewritten as 

(3.11) C{ft,K, t) = (/i - K)+ + e— dt'p{K, t'\ft) 
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Using the relation between the conditional probability t'\ft) and its gauge-transform P{s{f),t'\st), 
we obtain 



(3.12) C{fuK,r) = ((/* - K)+ + f P{s{K),t'\st)dt') 

Plugging the expression for P{s,t'\st) l|3.8() into (|3.12() and doing the integration over the time 
we obtain 0| 

(3.13) C{fuK.T) = e-^-{{h - K)+ + ^ 0„(s(i^))0„(sO^i^f^) 

A specific local volatility model will give an exact solution for a vanilla option if we can find the 
eigenvalues and eigenvectors for the corresponding Euclidean Schrodinger equation with a scalar 
potential. As examples of solvable potentials, one can cite the harmonic oscillator, Coulomb, 
Morse, Poschl- Teller I&II, Eckart and Manning-Rosen potentials. The classification of solvable 
scalar potentials was initiated by Natanzon )22| . This work provides the most general potential for 
which the Schrodinger equation can be reduced to either a hypergeometric or confluent equation. 
We will review in the following section the Natanzon potential, which depends on 6 parameters. 
We will later show that the Schrodinger equation can be doubled into a set of two independent 
Schrodinger equations with two different scalar potentials which transform into each other under 
a supersymmetric transformation. Moreover, if one scalar potential is solvable, the other one is. 
Applying this technique to the Natanzon potential, we will obtain a new class of solvable potentials 
corresponding to a new class of solvable diffusion processes. For these models, the solution to the 
KBS equation is given by a sum of hypergeometric functions. 

3.2. Solvable (super)potentials. In this section, we show that the Schrodinger equation can be 
formulated using supersymmetric techniques (see [H] for a nice review). In particular, the local 
volatility will be identified as a superpotential. Using this formalism, we show how to generate a 
hierarchy of solvable diffusion models starting from a known solvable superpotential, for example 
the hypergeometric or confluent hypergeometric Natanzon superpotential. 

3.2.1. Superpotential and local volatility. Let us write the Kolmogorov equation (|3.H|I in the fol- 
lowing way by introducing the flrst-order operator Ai and its formal adjoint A| ^ 

(3.1) 

with Ai = ds + W^^\s) , A[ = —da + W^-^\s). W^^^ is called the superpotential which satisfles 
the Riccatti equation 

(3.2) Q(i)(s) ^a.VFfi^s)- VF(i)'(s) 

Surprisingly, this equation is trivially solved for our specific expression for Q'^-* (|3.7|I (even with a 
drift M(/)!) 

■^In order to obtain the correct adjoint operator on R+, we impose the absorbing boundary condition <f)n{s = 
0) = 
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(3.3) 

^ ' ^ ' 2 ds V2a(i)(s) 

In particular, for zero drift, the local volatility function is directly related to the superpotential by 
(j(s) = e'^^ w(z)dz^ similar correspondence between the superpotential and driftless diffusion 
process has been found in ^| Moreover, if we have a family of solvable superpotentials 

^soivabie(^) ^ we Can always find an analytic solution to the Kolmogorov equation for any diffusion 
term a^^^ (s) by adjusting the drift with the relation (|3.3|) 

(3.4) m'^'(.s) = ^ - x/2.(i)(.)W^£L,,(.) 



Note that admits a zero eigenvalue if and only if the Kolmogorov equation admits a stationary 
distribution, 
distribution 



distribution. By observing that A\Ai(j)^Q '' = is equivalent to Ai^g^"* ~ 0, we obtain the stationary 



(3.5) <fl^\.s) = Ce-rw^''i^)d^ 

with C a normalization constant. Therefore, the stationary distribution will exist if the superpo- 
tential is normalisable (i.e. cp'j^^s) g L^). 

Example 3.2.1 (Coulomb superpotential and CEV process). The CEV process corresponding to 
cr(/) = o'o/'^ /i(/) = has the Coulomb superpotential W{s) = 2{i-f))s • 

Next, we define the Scholes-Black equation by intervening the operator Ai and A\ 

(3.6) El?^^^^\s) = 

(3.7) = (-a,2-Q(2)(s))4^)(s) 



We obtain a new Schrodinger equation with the partner potential 
(3.8) Q''^\s) ^ -dsW'^^^ - [W'^^^f 

Plugging our expression for the superpotential H3.3|l in (|3.8(l . we have 



In the same way as before, H2 admits a zero eigenvector (i.e. stationary distribution) if (j)Q (s) = 
Ce^ VK<^'(z)dz normalisable. 

Remark 3.2.2. In physics, the supersymmetry (SUSY) is said to be broken if at least one of the 
(1 2) 

eigenvectors (f>Q ' (s) exists. Otherwise, SUSY is said to be broken dynamically. 



10 



PIERRE HENRY-LABORDERE 



Now we want to show that if wc can solve the equation (|3.1() then we have automaticahy a solution 
to (|3.7|) and vice-versa. The SUSY-partner Hamiltonians Hi = A\^Ai and H2 = AiA\ obey 
the relations A\H2 = HiA[ and H2A1 = AiHi. As a consequence Hi and H2 are isospectral. 
More precisely, the strictly positive eigenvalues all coincide and the corresponding eigenvectors are 
related by the supercharge operators Ai and A\ : 

\> If Hi admits a zero eigenvalue (i.e. broken supersymmetry) , we have the relation 

(3.10) E^:^^=E^li; £;«==0; <t>^^\s) ^ Ce' r ^'^'i^)'^^ 

(3.11) 0i2)(.) = (4^_/i)-^Ai0i^ji(.) 

(3.12) 0W(s) = (i?^)-^40(^)(s) 

O If Hi (and H2) doesn't admit a zero eigenvalue (i.e. unbroken supersymmetry) 



(3.13) E^^)=Ell) 

(3.14) ,/.(2)(^)^(^(i))-i^^^(i)(^) 

(3.15) <j^^^\s)^{El)-'^AW{s) 

In the broken SUSY case, there are no zero modes and consequently the spectrum of Hi and H2 
are the same. One can then obtain the solution to the Scholes-Black (resp. Black-Scholes) equation 
if the eigenvalues/eigenvectors of the Black-Scholes (resp. Scholes-Black) are known. We clarify 
this correspondence by studying a specific example: the CEV process df = f^dW. In particular, 
we show that for /3 = |, the partner superpotential vanishes. It is therefore simpler to solve the 
Scholes-Black equation as Scholes-Black (rather than Black-Scholes) reduces to the heat kernel on 
. Applying a supersymmetric transformation on the Scholes-Black equation, wc can then derive 
the solution to the Black-Scholes equation. 

Example 3.2.3 (CEV with /3 = 2/3 and Bachelier process). We saw previously that the superpo- 
tential associated with the CEV process is given by 



(3.16) W'^^\s) ^ 



with the flat coordinate s = (iLq) ^i^d the potential H3.7|l 



2s(l-/3) 

(1-/3) 



(3.17) Q''H^)--.r. .,02 



4(1 -/3)2s2 



from which we deduce the partner potential (|3.9|l 



This partner potential corresponds to the potential of a CEV process df = f^dW with B given as 
a function of /3 by 
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B{B - 2) _ /3(2 - 3/3) 
(1-5)2 (1-/3)2 

Surprisingly, we observe that for /3 = |, Q'^'(s) cancels and the corresponding partner local 
volatility model is the Bachelier model df = dW for which the heat kernel is given by the normal 
distribution. The eigenvectors of the supersymmctric Hamiltonian partner H2 = —d^ to Hi are 
given by (with the absorbing boundary condition ^^(O) — 0) 



(3.20) ^^'^{s,E) = 



2)1 sin{\/Es) 



4:TTEi 



with a continuous spectrum M^. Applying the supersymmctric transformation (|3.12(l . we obtain 

2 
3 



the eigenvectors for the Hamiltonian Hi = —d^ + corresponding to the CEV process with (3 — - 



(3.21) cl)^'Hs,E) = E-^{-ds + -)(l)'-''\s,E) 

s 

1 , , rr: •. sin(\/Es) ^ 

(3.22) = ^ , {^VEcos{VEs) + ^- '-) 

VAttE'i s 

Plugging this expression in 13.13|l . we obtain the fair value for a vanilla option which can be 
integrated and written in terms of the cumulative distribution |19| 



e 



^-C{fu K, t) = {ft - K)+ + -^J^ '^^—^ ^'^''^ ^^"^^'^ ^) 



The fact that the CEV model with fi ~ depends on the normal cumulative distribution and is 
therefore related to the heat kernel on IR+ has been observed empirically by Here we have 

seen that it corresponds to the fact that the supersymmctric partner potential vanishes for this 
particular value of (3. 

3.3. Hierarchy of solvable diffusion processes. In the previous section we saw that the op- 
erators Ai and A| can be used to factorize the Hamiltonian Hi. These operators depend on the 
superpotential W^^^ which is determined once we know the first eigenvector (j)^^ [s) of Hi H3.5(l . 
We have assumed that Hi admits a zero eigenvalue. By shifting the energy E it is always pos- 
sible to achieve this condition. The partner Schrodinger equation (|3.7|l can then be recast into a 
Schrodinger equation with a zero eigenvalue 



(1) 



(3.1) H^.2)^ AiA\^ A\A2 + e[ 

where A2 = ds + W2{s) and A\ = ~ds + W2{s), 



(3.2) M^(2)(,)-1^^-^^H^) 



2 ds V2ct(2)(s) 
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We have introduced the notation that in En , n denotes the energy level and (m) refers to the 
m*^ Hamiltonian By construction, this new Hamiltonian H2 = ^2^2 + e[^'' is solvable as 

AiA\ is and the associated diffusion process with volatility cr^^-* and drift /i^^^ satisfy (|3.2|l . The 
superpotential W'^'^\s) is determined by the first eigenvector of H2, 4'^\s), 



(3.3) 

We can then apply a supersymmetric transformation on H2. The new Hamiltonian H3 can be 
refactorised exactly in the same way we did for H2. Finally, it is not difficult to see that if Hi admits 
p discrete (normalisable) eigenvectors, then one can generate a family of solvable Hamiltonians Hm 
(with a zero-eigenvalue by construction) 



(3.4) H„, = AlA„, + eI^J_, = -dl + g,„(s) 

where A,„ = dg + Wm{s). This corresponds to the solvable diffusion process with a drift and a 
volatility such that 



(3.5) VF„(.)= '^^-^"^l^'-^-H^) 



ds 2 ds V2cr(™)(s) 

The eigenvalues/eigenvectors of Hm are related to those of Hi by 



(3.6) i?i")=i?iV) = ... = i?«_, 

(3.7) 01™) = {Ei'l^^, - eIZ2)-^ ■ ■ ■ (eUL-i 4")-^ A,„_i . . . 

In particular, the superpotential of _ff,„ is determined by the (m— 1)"* eigenvector of -ffi, 0[^Li(s) 



,0 w ( \ dlniAm .^i ■ ■ ■ Ai^b'mliis)) 
(3.8) Wm(s) = 



(3.9) 



ds 



2 ds V2a(™)(s) 



Consequently, if we know all the m discrete eigenvalues and eigenvectors of Hi, we immediately 
know all the energy eigenvalues and eigenfunctions of the hierarchy of m — 1 Hamiltonians. In 
the following we apply this procedure starting from a known solvable superpotential, the Natanzon 
superpotential. 



3.4. Natanzon (super)potentials. The Natanzon potential is the most general potential 
which allows us to reduce the Schrodinger equation H3.6() to a Gauss or confluent hypergeometric 
equation (GHE or CHE). 
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3.4.1. Gauss hypergeometric potential. The potential is given by 

with i?(z) = r2Z^ + riz + ro > and S'(z) — S2Z^ + siz + so (two second order polynomials). 
The z coordinate, lying in the interval [0, 1], is defined implieitly in terms of s by the differential 
equation 

(3 2) = 2z(l-z) 

Example 3.4.1. The hypergeometric Natanzon potential includes as special cases the Posch- Teller 
potential II Q(s) = A^Bsech{^^)'^ + Ccsch{^^)'^ for ro = r2 = and the Rosen-Morse potential 

Q{s) = A + Btanh{^) + Csech{^Y, ri = r2 = with A and B constant. 

By construction, the solution to the Schrodinger equation with a GHE potential is given in terms 
of the Gauss hypergeometric function F(a,(3,"f, z) 

(3.3) E) = {z'r^zi (1 - z)"^^^F(a, /3, 7, z) 

where F{a, (3,j, z) satisfies the differential equation P 



(PE , , , ,dE 

■^ + (7-(a + /3+l)z)- 

and the most general solution to this equation (|3.4fl is generated by a two-dimensional vector space 



(3.4) z{l - z)— + (7 - (a + /3 + l)z)— ~a[3z^Q 



(3.5) F{a, /?, 7, z) = C12 Fi (a, /?, 7, z) + 022^-^2 Fi (a - 7 + 1, /3 - 7 -f 1, 2 - 7, z) 

with 2 -Fi(a, /3, 7, z) satisfying 2 Fi (a, /3, 7, 0) = 1 and ci and C2 two arbitrary coefheients. The 
parameters a, /?, 7 depend explicitly on the eigenvalue E by 



(3.6) l-{a- (3f ^r2E + S2 

(3.7) 27(a-h^-l)-4a/3 = ri£; + si 

(3.8) 7(2 ~ 7) = roF + so 

In fact, one can show that E satisfies a fourth-order polynomial and find E as an explicit function 
of a, /?, 7 jJT]. By imposing the condition that the eigenvectors are normalisable (i.e. belong to 
i^([0, 1])) we obtain the discrete spectrum En and can determine the coefficients ci and C2. We 
impose conditions on a, /3, 7, ci and C2 such that 

(3.9) dz:^z-/-2(i _ z)-^+«+'3-iF2(a, /3, 7, 2) < 00 
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Looking at the asymptotic behavior of 2^"i(a, /?, 7, near z — and z = 1 ^ ^, wc obtain the 
following conditions 





ci = 0, 7 < 1 or C2 = 0, 7 > 1 


ro = 0, ri ^ 


ci = 0, 7 < 2 or C2 = 0, 7 > 


?'o = 0; ri = 


ci = 0, 7 < 3 or C2 = 0, 7 > — 1 



Table 2. condition z ~0 



ci = 


a-leN*,a + /3-7<0or-l-Q; + 7eN*,a + /3-7>0 


C2 = 


-aeW,a + f3-j>0oia--feW , a + f3 - j < 



Table 3. condition z = 1 



We have the discrete eigenvalues (a„ = —n ; n € N*) 

(3.10) 2n + 1 = -v/l - roEn -so + y/l- r2S„ - S2 - Vl - (''o + + r2)-B„ - {sq + si + S2) 

(3.11) ^„(s) = (z')"^2*(^- l)^^"^^*^^i^(an,/3«,7n,^) 

(3.12) 7„ = 1 + v/l - roi;n - So 

(3.13) a„ - /3„ = -V^l - ''2£^n - S2 

(3.14) a„ + /3„ = 7„ + ^1 - (ro + ri + r2)£'„ - (so + si + S2) 
and the (normalised) eigenvectors 



(3.15) Ms) = i3„(^')"^^*(l - ^) ''"'"^ - 2z) 

with i?„ = [((^) + (^) - (^))£M!H^W^]-i Pi--^-+'^--)(2z - 1) the Jacobi 
polynomials. 

3.4.2. Confluent hypergeometric potential. A similar construction can be achieved for the class of 
scaled confluent hypergeometric functions. The potential is given by 

(3-16) Q{s) = , (—-7 ^77^ ^2) 



R{z) 4 R{z) ' R{zY 

with R{z) = r22;^ + riz + ro > and S{z) = S2Z^ + siz + so. The z coordinate, lying in the 
interval [0,cx)[, is defined implicitly in terms of s by the differential equation 

dz{s) 2z 



(3.17) 



ds ^/R{z) 



Example 3.4.2. The confluent Natanzon potential reduces to the Morse potential Q{s) = 



-l+SQ+sie n/^ +S2e 



for ri = r2 = 0, to the 3D oscillator Q{s) = — ^5—^ + 77 + ^'^'^ ro = r2 = and to the Coulomb 
potential Q(r) = -'''"-^l^^:'-''''' for ro = n = 0. 

.) r(7) ({-1 + z)—^+~< + 
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By construction, the solution to the CHE potential is given in terms of the confluent hypcrgeomctric 
function F(a, /?, 7, z) 

(3.18) ^{s,E) ^ z{s)ie-^{z'{s))-^F{a,P,"f,ujz{s)) 
The parameters u>, 7, a depend explicitly on the eigenvalue E by 

(3.19) = -ra-B - S2 

(3.20) 2w(7-2a) =ri£ + si 

(3.21) 7(2-7) =roS + so 

Note that 0(z) = F{a, (3,j,ljz{s)) satisfies the differential equation P 



(3.22) z(j)"{z) + (7 - wz)0'(z) - uja(l){z) = 

and the most general solution to (|3.22|) is generated by a two-dimensional vector space 

(3.23) F{a, 7, uz) = aM (a, 7, wz) + C2(wz)^"'^A/(l + a - 7, 2 - 7, wz) 

with 7\/(a,7,a;z) the M-Whittaker function and ci and C2 two arbitrary coefficients. By imposing 
that the eigenvectors are normalisable ((i.e. belong to i^(R+))) , we obtain the following conditions 
(see Table 2 &5) which give the spectrum E and the coefficients ci and C2 ^ 



a > 2 


no condition 


a < 2 


ci 0, -1 -q: + 7 e N* 
or C2 = 0, -a e N* 



Table 4. condition z = 00 



We have the discrete eigenvalues Q!„ = — ?i ; n e N 



(3.24) 7„ = 1 + v/1 - roi;„ - so 

(3.25) [J„ = \/-r2En - S2 

(3.26) 2n + l = , ^ - y/l - rpEn - sp 
and the (normalized) eigenvectors 



(3.27) ^„(s) = ^z{sy-^e-^{z'{s))-'^LZ--\oj^zis)) 

with LJ^"~^(z) the (generalized) Laguerre polynomial. In the following, we have listed classical 
solvable superpotentials and the corresponding solvable local volatility models and solvable instan- 
taneous short-rate models (Table 5 & 6). 
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Superpotential 


W{s) 


£(£) 

T-IJ 


Shifted oscillator 


as + b 




Coulomb 


« + ! 




Morse 


a + be-'^' 




Eckart 


acoth{as) + b 




Rosen-Morse 


ath{as) + 6 


g2&s+^/n(cosh(Q)) 


3D oscillator 


as + \ 


gas^+2Wn(s) 


P-T I a > 26 


atan{as) + bcotg{as) 




P-T II a > 26 


ath{as) + bcoth{as) 


^^alog(co3h(Qs)) 1 blog(sl„h(as))^ 



Table 5 . Example of solvable superpotentials-Local Volatility 



one-faetor short-rate model 


SDE 


Superpotential 


Vasicck-HuU- White 


dr = k{e- r)dt + adW 


Shifted Osc. (a = f , 6 = 


CIR 


dr = k{e ~ r)dt -f a^/rdW 


3D Osc. (a^|,6=i-f|^') 


Doleans 


dr — krdt + ardW 


(Constant {W{s) - -'^''^''^ ) 


EV-BK 


dr = r{ri — aln{r))dt + ardW 


Shifted Ose. (a - ^, 6 - v^2(-2')+^^) ) 



Table 6. Example of solvable one-faetor short-rate models 



3.4.3. Natanzon hierachy and new solvable processes. We know that the Natanzon superpotential 
is related to the zero-eigenvector 



(3.28) Wnat = -dMMs)) 

and the corresponding supercharge A is 



2z(l-z) _ 70 _ (l + ao + /3o-7o) ao/3aF{l + a^,! + pQ,! + jo, z) , z'-lz"{s). 



VR 2z 2(z-1) 7oF(ao,/3o,7o,2) 2 

With a zero drift, the Natanzon superpotential corresponds to the diffusion process H3.3|l 



(3.29) ^(i)(,) = ^L 

with (Tq^'' a constant of integration. Applying the results of section (3.2), we obtain that the driftless 
diffusion processes 



(3.30) 



'(.) 



Am) 



A, 
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are solvable ()3.9|l . Using the fact that the (discrete) eigenvector 4'lJ_i{s) is a hypergeometric 
function and that the derivative of a hypergeometric function is a new hypergeometric function ^ 
the action of A^-i ■ ■ ■ Ai on 0mi (s) will result in a sum of (m — 1) hypergeometric functions, thus 
generalizing the solution found in [2]. 



4. Gauge-free stochastic volatility models 

In this section, we try to identify the class of time-homogeneous stochastic volatility models which 
leads to an exact solution to the KBS equation. We assume that the forward / and the volatility 
o are driven by two correlated Brownian motions in the risk-neutral measure 



(4.1) dft ^ atC{h)dWi 

(4.2) dat = b{atJt)dt + a{atJt)dW2 

(4.3) dWidW2 = pdt 

with initial condition oq — a and /o = /• 

Using the definition for the connection H2.11|) . we obtain the Abclian connection 

(4.7) A ^ (^(^ - lojnC-)) 
Then, the field strength is 

(4.8) J^af = daAf-dfAa 



We will now assume that the connection is flat, Taf = 0, meaning that the connection can be 
eliminated modulo a gauge transformation. In this case, the stochastic volatility model can thus 
be called a gauge-free model. This condition is satisfied for every correlation p if and only if 

(4.10) daMa)^df^ 

(4.11) da -dl- + —d}- = 



^2Fi(a,l3,-,,z) = ^2Fi(a + l,l3 + l,-/ + l,z) and M'(a, /3, z) = !lM{a + 1, /3 + 1, z) 
7 

(4.4) = df — 

(4.5) -^"-^(-T^^^^) 
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Moreover, if we assume that a°'{a) is only a function of a (this hypothesis is equivalent to assuming 
that the metric admits a Killing vector), the model is gauge- free if and only if 

(4.12) '- = -.d/- + aHf)~^-^dnf) 



with 0(/) satisfying 
(4.13) df^if) 



a 2 a 2 ' ^ ' J cr(a') 

df{Cd)C) fa'da' 



a{a') 



This last equation is equivalent to C{f)djC{f) = (3 with /3 a constant and = 7 a constant 
function. For /3 = 0, the above equations simplify and wc have 



(4.14) C(/) = /i/ + 

(4.15) 6(a) = a(T(a)(7 + iaa^) 

2 a 

with ^, 1/ , 7 constant. The gauge-free condition has therefore imposed the functional form of the 
drift term. When the volatility function is fixed respectively to a constant (Heston model), a linear 
function (geometric Brownian model) and a quadratic function (3/2-model) in the volatility, one 
obtains the correct (mean-reversion) drift ® (see Table 7) 



name 


cr(a) 


SDE 


Hcstoii 


a{a) = -q 


dv = ^/5{2v■^ + ri{ri - l))dt + 2i'iVSy/^dW2 


Geometric Brownian 


cr(a) = 'qa 


dv = VS{2r]^vi + ifv)dt + 2\/5rivdW2 


3/2-modcl 


a {a) = rja^ 


dv = 2^/5r]{r] -{■ j)v'^dt + 2VSr]v^dW2 



Table 7. Examples of Gauge free stochastic volatility models with df = S{fj, + I'f jy^dWl. 
The gauge transformation eliminating the connection is then 

(4.16) A(/, a) = -P^lniCif)) + apf^ ^ + (/^a 
Finally, plugging the expression for C{f) and h{a) into (|2.12l) 

(4.17) Aa^ + B(j^da - 

a 

The Black-Scholes equation for a Vanilla option C{a, f,T = T — t) (with strike K and maturity T) 
satisfied by the gauge transformed function C'{a, /, r) = e-'^^^'°'''C{a, /, r) is 

^In order to obtain the correct number of parameters, one needs to apply a change of local time dt = 5dt' , 
dWi,2 = VSdW{2 

and B arc two constants given by 

A = ^{-^p"dfC + ap'^ff + ^pffC^d]ln{C) + p{^^pffdfC + ap''f)iapff - ^dfC) + i(p^^a - ^dfCf 
1 .,dfC. 



, wc find that the function Q is ^ 
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(4.18) drC'ia, /, r) = AC (a, /, r) + Q{a)C'ia, /, r) 

with the initial condition C'{aJ,T = 0) = e^^f'"-\f - K)+ . In the coordinates q{f) = 
and a, the Laplace-Beltrami operator is given by 

(4.19) A^^{-d^ + 2pda,+da-da) 

la a 
Applying a Fourier transformation according to q, C'{q, a, r) = TC {f, a, t), we obtain 



(4.20) drC'{a,q,T) = —{-k''- + 2tkpda + da-da)C'{a,q,T) + Q{a)C'{a,q,T) 

la a 

with the initial condition C'{a,q,T = 0) = J^[e^'^^'°''' {f — K)^]. Using a spectral decomposition 
C'{a,q,T) = J2n'l^rik{a){C'{a,q,T = 0),(j)nk)e~^"'''^ (with (., .) the scalar product on L^), the 
eigenvectors 0„fc(a) satisfy the equation 



(4.21) - EnkKkia) = -7r{~^^- + '^ikpda + da-da)(f)nk{a) + Q(a)0nfc(a) 

la a 

This equation H4.2H) can be further simplified by applying a Liouville transformation consisting in 
a gauge transformation and a change of variable [211 

(4.22) Vnfe(s) = (-)^e'''''^°^'"V„fc(a) 

a 

(4.23) ^ 



da '^{o-) 
4'nk{s) satisfies a Schrodinger equation 

(4.24) V"fc(s) + {E^k - J(s))V„fc(.s) = 

with the scalar potential 



fcV 4 a4 fc2 ^2 _ 3 ^ ^2 ^,(^^Y ^ 2 Q cr(Q) ((j'(q) - a a" [a)) 1 

^ ^ ^ 2^ 



and where the curly bracket denotes the Schwarzian derivative of a with respect to s 

The 2d partial differential equation corresponding to our original KBS equation for our stochastic 
volatility model has therefore been reduced via a change of coordinates and gauge transformations 
to a Schrodinger equation with a scalar potential J(s). The stochastic volatility model is therefore 
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solvable in terms of hypergeometric functions if the potential J(s) belongs to the Natanzon class. 
The solution is given (in terms of the eigenvectors ipnk) by 

(4.26) C(a,/,r)=e-^('^'^).F-M^^nfc(s(a))(.F[e^(^'''^)(/-if)+],^„fc(s(a)))e-^-1 

n 

Let us examine classical examples of solvable stochastic volatility model and show that the corre- 
sponding potentials J(s) correspond to the Natanzon class (see Table 8). 



name 


potential 


J(,s) 


Heston 


3D osc. 


Jis) = 


-3+4 B s-' ?7+s* (2 A+k'' 








Geometric Brownian 


Morse 




-7f+4e-^^-^ (2A+fe2(- 
8 




3/2-model 


Coulonibian 


J{s) 


8 A+j;(-8B+J))+4fe^ (- 

4 ^2 jj2 


+ 



Table 8. Stochastic volatility models and potential J(s) 



We present here a new example of solvable stochastic volatility model which corresponds to the 
Posh- Teller I potential. 



Example 4.0.3 (Posh- Teller I). For a volatility function given by (j{a) = a + ija^, the potential is 
given by a Posh- Teller I potential 




5. Conclusion 

We have shown how to use supersymmetric methods to generate new solutions to the Kol- 
mogorov & Black-Scholes equation (KBS) for one-dimensional diffusion processes. In particular, 
by applying a supersymmetric transformation on the Natanzon potential, we have generated a 
hierarchy of new solvable processes. Then, we have classified the stochastic volatility models which 
admit a flat Abelian connection (with one Killing vector). The two-dimensional KBS equation 
has been converted into a Schrodinger equation with a scalar potential. The models for which the 
scalar potential belongs to the Natanzon class are solvable in terms of hypergeometric functions. 
This is the case for the Heston model, the geometric brownian model and the 3/2-model. A new 
solution with a volatility of the volatility a{a) = a + rjo^, corresponding to the Posh- Teller I, has 
been presented. 
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